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We present a continuum phase-field model of crack propagation. It includes a phase-field that is 
proportional to the mass density and a displacement field that is governed by linear elastic theory. 
Generic macroscopic crack growth laws emerge naturally from this model. In contrast to classical 
continuum fracture mechanics simulations, our model avoids numerical front tracking. The added 
phase-field smoothes the sharp interface, enabling us to use equations of motion for the material 
(grounded in basic physical principles) rather than for the interface (which often are deduced from 
complicated theories or empirical observations). The interface dynamics thus emerges naturally. In 
this paper, we look at stationary solutions of the model, mode I fracture, and also discuss numerical 
issues. We find that the Griffith's threshold underestimates the critical value at which our system 
fractures due to long wavelength modes excited by the fracture process. 

PACS numbers: 62.20.Mk, 46.50.-|-a 



I. INTRODUCTION 

The study of fracture is usually approached using 
mathematical descriptions and numerical simulations 
based on empirical observations. Finite element methods 
are commonly used to investigate the behavior of frac- 
tured materials on a large scale, where the crack growth 
lawsjl} H (that is, velocity and direction of the crack for a 
given stress field near the tip) are introduced empirically. 

We present a continuum description starting from ba- 
sic theoretical assumptions. We introduce a phase-field 
model, originally used to describe thermodynamic phase 
transitions and widely used to model solidification and 
combine it with a displacement field. In contrast to other 
phase-field models of fracture ^ and interfacial motion 
in the presence ofstrain|,|l, our phase field is conserved, 
representing the density of the material. Bhate et al. ^ 
study a conserved order-parameter phase-field model in 
the context of stress voiding in electromigration; their 
dynamics is rather different from ours, since their elastic 
deformations are quasi-statically relaxed (the limit in our 
theory of 77 ^ 0, see below). 

The phase-field serves two main purposes. First, it 
smears out any sharp interfaces, facilitating numerical 
convergence. Second, the model gives equations of mo- 
tion for the material rather than the boundaries, thus 
we avoid dealing with a moving boundary value prob- 
lem which would require numerical front tracking. One 
of our main goals is to find macroscopic fracture laws. 
In our model, these laws emerge naturally from the dy- 
namics of the fields. See Fig. |l| for a three dimensional 
representation of the phase-field in a fracturing sample. 

One incentive to use a conserved-order parameter is 
simply that density is conserved microscopically (apart 
from applications where etching or sublimation is impor- 
tant). In general, a non-conserved phase-field will give 
a non-zero velocity even for a straight material interface 




FIG. 1: A three dimensional representation of the phase-field 
in a fracturing sample. The vertical axis shows the value 
of the phase-field 4'{x,y), where = 1 is unstrained material 
and = is vacuum. This example corresponds to the fourth 
contour in Fig. M. The values of x and y are given in units of 
w/h (see App. 



y. This could be remedied by tuning the free energy 
so that the material and vacuum have the same energy 
density, but then the strained region around the crack tip 
would evaporate. Conserving the phase-field also gave us 
insight into how to properly implement the conservation 
laws (see Sec. |l|). Another option would have been to 
add a non-conserved field, such as damage or dislocation 
density, in addition to our conserved mass. This would 
add complexity without ameliorating the numerical chal- 
lenges presented by the conservation law. In future work 
we intend to introduce such non-conserved state variables 
to model plastic fiow. 

The next section gives an outline of the theoretical 
model, presenting the main equations. We then inves- 
tigate some of the stationary solutions analytically, and 
discuss their consequences. This is followed by a brief 
presentation of the numerical implementation. We then 
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measure the crack growth velocity as a function of ex- 
ternal stress and explore the fracture threshold of our 
model. We conclude with suggestions for future work. 



II. THE FRACTURE MODEL 

The model consists of a phase-field cj) and a displace- 
ment field u. The former field is interpreted as the nor- 
malized mass density, and typically has values between 
zero and one. The latter field, through its derivatives, 
represents strain in the material. The model is based on 
a free energy T . The equations of motion locally conserve 
the density cf), moves it under the flow-field u, and evolves 
both (f) and u in the direction of the net force from the 
free energy: in particular they are constructed so that 
dT /dt < 0. The free energy is given by the integral 
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where 



(1) 



(2) 



The first term in Eq. (|^) is a gradient term, energetically 
penalizing spatial fluctuations in the phase-field. The 
first term in Eq. (^) is a Ginzburg-Landau double well po- 
tential, favoring values of at zero and (f>s [e] = 1 — emm 
(using the Einstein summing convention), representing 
the two phases vacuum and solid, respectively. If the ma- 
terial is completely unstrained the solid value is 4>s [e] = 1 , 
otherwise this value is either higher (for a compressed 
material) or lower (for a stretched material), where emm 
is the the density change for small strain. The factor 
4>s[e\ — 4> can be thought of as a density of vacancies 
or interstitials. The parameter h controls the height of 
the energy barrier between the vacuum and solid phases. 
The ratio of w and h controls the width of the solid- 
vacuum interface, that is the width of the transition from 
= 05 [e] to = 0. 

The next term is the elastic strain energy density £ [e] . 
The elastic energy is calculated from the strain tensor e, 
and is given by 



£[e] 



(3) 



For a homogeneous, isotropic material, the tensor Cijki 
can be described by the two Lame constants jjL and A 
through 



kl 



(4) 



where /i is the shear modulus and A is proportional to 
the bulk modulus. This gives 



(5) 



In two dimensions, we get 



The strain field eij is related to the displacement field by 
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Note that in this definition of the strain field, we are ig- 
noring geometric nonlinearities|^, which are important 
for large rotations. According to Eq. (|^), the divergence 
of the displacement field is just the trace of the strain, 
V • u = Emm- The displacement field u(x) is defined 
in the deformed or Eulerian coordinate system, which 
means that x describes a location in space. (In the unde- 
formed or Lagrangian description, x would describe the 
location of the material before the displacement is taken 
into account; Lagrangian coordinates are usually used in 
finite element calculations.) The Lame constants are con- 
nected through the Poisson ratio 1/ by A = 2/Ltj// (1 — 2v), 
see Appendix In the case of plane strain, the addition 
of the V • u term in the double well potential turns out 
to be crucial to preserve this relation. Since the elastic 
energy f [e] is only defined in the material (that is, where 
7^ 0), the elastic term is multiplied by a factor of 0^; 
thus the strain energy will go to zero in the vacuum. 

The equations of motion we have chosen for the phase- 
field (j) and displacement field u are overdamped and Eu- 
lerian, moving the fields along the direction of net force. 
The time derivative is thus proportional to the force on 
the field. Physically, our model might describe fracture of 
a colloidal crystal, or "atoms in molasses" . We are there- 
fore intermediate between quasi-static fracture (where 
the crack evolution is calculated from the static strain 
of the current configuration) and dynamic fracture (with 
inertial effects and wave reflection at the boundaries). 
Specifically, our equations of motion are 



dt 



= -V - J 
du 



J = -DV 



5T 



dvL 



(8a) 
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(8b) 



where 77 and D are the viscosity and the diffusion con- 
stant, respectively. Note that Eq. ( ^a|) is the continuity 
equation. This means that total 0, or mass, is conserved. 
The first term in J is a diffusion term, while the sec- 
ond term makes sure that the mass follows the motion of 
the displacement field. The total variational derivative 
DjF/ISu in Eq. (^) can be found by first noting that 
a small change Au in u results in a change in 0. The 
new value of at a point changes due to two effects: 
(1) a gradient in dragged by a distance Au changes it 
by — (V0) • (Au), and (2) a divergence in Au causes a 
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change in density — 0V • (Au). Together these combine 
into the net change A(/)[Au] = —V • {(pAu). This is the 
continuity equation, where 0Au is the flux of (j>. The 
total change in the free energy is then 
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.Au) 



dV 



— • Au - V 
ou 

^ • Au + (/)V^ • Au ) 
ou ocp 



(9) 



We assume that the boundary terms vanish in the inte- 
gration by parts. 

With the equations of motion (||), the evolution of the 
fields are overdamped and the free energy is decreasing 
in time, 
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(10) 



To solve for the fields, the functional derivatives need 
to be calculated explicitly. They can be written in the 
convenient form 



ST 

Su. 
ST 



dg 



d(t> 



(11a) 



(lib) 



Finally, it should be pointed out that all the equations 
can be made unitless by rescaling all the quantities in- 
volved. For the present model it is convenient to use 
w/h as the unit length, as the unit energy density, 
and l/jy as the unit diffusivity. This corresponds to set- 
ting w = 1, /i = 1 and ry = 1 in Eqs. (jT]), (js^), and ([u]). 
See Appendix ^ for information on reduced units for the 
quantities used in this paper. 



that V\STIS4)\ = and ST/Su = 0. Physically ST/Stj) 
can be considered a chemical potential, and a non-zero 
gradient will result in a flow of material. This means 
that the chemical potential p is a constant. Thus to flnd 
a stationary solution, we have to require that 



ST 



P 



and 



ST 
Su 



0. 



(12) 



We will find the stationary solution of a single straight 
interface running perpendicular to the a;-direction be- 
tween the solid and the vacuum. We will therefore as- 
sume that (j) and e^x only vary with respect to x, that eyy 
is constant, and that exy — £yx = 0. Combining Eqs. (|l|) 
and (|l^) with the requirement in Eq. (p^), we get 



dg 
dg 



Integrating Eq. (13b) gives 

dg 

dCxx 



vv 



C. 



(13a) 
(13b) 
(13c) 

(14) 



For non-zero C, the vacuum phase {((> ~ 0) can be shown 
to be unstable, so to get an interface we must set C = 0. 
Solving Eq. dl4) with respect to e^x gives 



exx = (1 - A) [1 - ,^ - (1 + 2\)eyy] 



where 



A- 



A + 2/i+ 1/2 



(15) 



(16) 



Next we multiply Eq. ( p^ with dx4>. Using Eq. (|lj) 
and reme mbering that dx^yy = 0, we can then rewrite 
Eq. (|l3a| ) as 



dx 



0. 



Upon integration, this becomes 
1 



(17) 



(18) 



where 



III. STATIONARY SOLUTION 

This section will calculate the profile of a stationary 
straight interface between the solid and vacuum phases. 
Before delving into the details, consider Eqs. (^). A sta- 
tionary solution means that <j) = and u = 0. Unless we 
have the trivial solution where (j) is constant, this implies 



T[4>]=g[c^,e[<l>]]-p^-q, (19) 

q is the constant of integration, and Eq. (|l^) has been 
used to write g[(j), e] as a function of (j) only {cyy is consid- 
ered a constant). Here T[(j)] can have two minima, giving 
the solid and vacuum densities. 

Notice that Eq. (|8|) looks like the Hamiltonian of a 
classical particle at position and time x, where the 
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first term is the kinetic energy and —T[(j3\ is the potential 
energy. If we want a solution that starts at a small and 
constant (j) at x = —oo, and ends at a larger constant (j> at 
X = oo, then T[0] must have two stationary points with 
respect to </>, and these two points must have the same 
value. Since T[0] is fourth order in 0, it can be written 
in the general formp7| 



2? 



(20) 



Comparing this to Eqs. (||) and ([Tg|), we find that|T6 
B"^ = A/2 and 



»1,2 



where 



{A 



2k{2A-k) 2 
1-A VV 



2A 



K = A(1 + 2A)-2A. 
To second order in Cyy this gives 



k{2A-k) 



2A2(1 



A) 

k{2A- k) 



^ A^^^ 2A^i-Afyy 



(21) 
(22) 

(23a) 
(23b) 



Note that eyy = gives (pi — and 02 = 1- 

Inserting Eq. (|2Q[ ) into Eq. (|l^) and solving for <j) gives 

= ^ tanh [5B{x - c)] + (/-o (24) 

where 
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(25) 



and the constant of integration c determines the location 
of the interface. In Fig. ^ we plot e^x and as a function 
of position according to Eqs. (|l5|) and (p^) for A = 2/i = 



2, c = and eyy = 0. Notice that the form of Eq. (|2i 
implies 




0.5 - 



-5 5 

X (in units of w/h) 



FIG. 2: A plot showing the stationary solution when eyy = 0, 
with X — 2fi = 2 and c = 0. Note that the strain e^xix) stays 
non-zero as (^(a;) — ^ 0; this is acceptable, since the strain does 
not contribute to the free energy in this limit (this is not true 



when eyy ^ 0) . 



The approximation is valid if the interface is far away 
from the boundaries of the sample, which means that 
Xi,2 > 1/5B. 

As a prelude to the section on numerical implemen- 
tation, we should point out some shortcomings of the 
current definition of the free energy. As seen above, the 
limiting value of (f) in the vacuum side of the interface is 
not zero if eyy ^ 0. Strictly speaking, there is no longer 
a vacuum, but a gas filling the voids of the cracks. What 
makes this troublesome is that this "gas" can support 
shear forces. We have therefore chosen to do the nu- 
merical simulations using eyy — when measuring the 
fracture threshold. To improve the model for more com- 
plex runs, the free energy could be adjusted to assure 
that the value of (f> in the limit of no material is zero for 
all stationary solutions. 



IV. NUMERICS 



A. Implementation 



P 



2B^ 



and 



B' 



q 



(26) 



Consider a sample that is X wide and y tall, with 
an interface perpendicular to the x-direction. Let X = 
Xi + X2, where Xi is the amount of the sample that has 
4> < (t>o (vacuum), and X2 corresponds to > 0o (solid). 
Then the free energy per unit length is 



1 



2T[4>] 



q dx 



X 



U^B+p{(l)iXi+^2X2) + qX . 
o 



(27) 



We have implemented Eqs. ( pi] ) for a plane strain 
system. Thus we can perform our simulations on a 
two-dimensional uniform structured grid with periodic 
boundary conditions in both directions. The periodic 
boundary conditions allows the use of Fourier methods. 
To increase stability, we implemented a semi-implicit 
scheme pO|. The linear terms can be solved analyti- 
cally in Fourier space, which increases efficiency con- 
siderably. Specifically, at each timestep we first inte- 
grate the nonlinear terms using an explicit Euler scheme 
before multiplying with the factor exp(— dtV^), where 
is the discrete version of the Laplace operator in 
k-space. In our case, this operator is equal to = 
X)i=i 2{[2 cos{kiAxi) — 1]/ (Axi)"^}. The exponential fac- 
tor represents the analytical solution to the linear part 
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FIG. 3: A strained material with a double-ended crack. The 
hollow arrows indicate the loading direction and the dashed 
lines are the periodic boundaries. 



of the time derivative, (p — —V 0. 

In Fig. H we show the setup for a double-ended crack 
under mode I loading that we use in our numerical sim- 
ulations. (See Fig. |l| for a three-dimensional representa- 
tion.) 

The system is initially strained in the x-direction (the 
horizontal direction in Fig. ^ with a uniform constant 
strain e^x- Numerically, the strain is represented through 
the spatial derivatives of the displacement field u, which 
means that there is an inherent discontinuity in the strain 
field at the left and right boundary in Fig. ^. This prob- 
lem has been resolved using "skew-periodic" boundary 
conditions. In essence, we identify u on the left with 
u + Au on the right. 

The initial phase-field is set to the constant value that 
minimizes the free energy for the uniform initial strain. 
A circular hole is inserted in the middle by removing 
mass (that is, by tapering the phase-field to zero) in a 
circular area in the center. A crack will grow if the strain 
exceeds the fracture threshold. We will later compare 
this threshold to the Griffith's criterion. 

The fourth-order gradients in our evolution equations 
(^ make this problem numerically challenging: most 
simple algorithms will go unstable at a time step which 
goes as the fourth power of the grid spacing. Roughly 
speaking, the time step must remain smaller than the 
time it takes information to pass across the footprint 
of the numerical stencil (the size of the region used to 
calculate the gradient terms). To make for an efficient 
algorithm, we paid careful attention to minimizing this 
footprint area: in so doing, we found that it was impor- 
tant to pay close attention to the locations of the various 
terms with respect to the numerical grid. (For example, 
the asymmetric forward derivative of the phase field is 
located at the midpoint of the bond between two grid 
points). Ref. contains further details and sugges- 
tions. Notice that the footprint was also reduced through 
our use of the semi-implicit scheme mentioned earlier, 
leaving only third-order gradients to be solved numeri- 
cally. 



B. Results 

Consider a block of material with a small flaw in it. If 
an increasing strain is applied, then at a certain point the 
material will fail, and a crack will form. Griffith sug- 
gested that this would happen when the strain energy 
released is just enough to form the two newly created 
crack surfaces. Below we explore how our numerical sim- 
ulations compare with Griffith's simple model. 

In most materials, the actual fracture threshold (that 
is, the strain energy at which the material fractures) is 
higher than the Griffith's threshold. The reason is that 
the strain energy is converted not only to surface energy, 
but also to plastic deformation, sound emission, and heat. 
Our model has no plastic deformation, but we will see 
that some of the strain energy is lost through long wave- 
length emission similar to that of phonons in dynamic 
fracture. 

When measuring the fracture threshold experimen- 
tally, the load or displacement on a specimen is mono- 
tonically increased until it breaks. This is difficult to do 
in our simulations, so instead we run a series of simu- 
lations, each with a different strain as initial condition. 
We then measure the crack tip velocity as a function of 
strain energy per unit length stored in front of the crack, 
and define the fracture threshold as the energy where the 
velocity extrapolates to zero. 

To measure the crack tip velocity we simulated a dou- 
ble ended crack (see Fig. |) on a (A", 3^)=(200,200) grid 
(Ax — Ay = 1) with periodic boundary conditions. The 
Lame constants were chosen as A = 2/i = 2, which corre- 
sponds to v — 1/3. Initially, each simulation was given a 
uniform strain e^x in the x-direction and no strain eyy = 
in the y-direction. The phase-field 4> was uniformly set 
to the value that minimizes the free energy integrand in 
Eq. (pi) for a constant strain|l9l: 
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-0.[e] + -y^![eF32£R, 



(28) 



To initiate the crack, a circular hole of radius 10 was in- 
serted in the center of the sample; here, a "hole" means 
that we set to zero, being careful to make the edges 
smooth in order to avoid numerical instabilities. After 
an initial transient period, the crack would start to grow 
in the y-direction at a uniform rate until it reached the 
boundaries. Here, the crack would sense its periodic im- 
age and speed up as the two crack tips coalesced. Fig. ^ 
shows part of the grid with the (p — 0.5 contours as the 
double ended crack grows in the y-direction. 

To measure the crack tip velocity, we needed to track 
the crack tip. It turns out that the free energy density has 
a peak in the tip area, so we decided to define the location 
of this peak as the crack tip position. The velocity could 
then easily be found by finding the slope of the curve 
describing the crack tip position as a function of time 
(we used the tip growing in the positive y-direction in 
Fig. ^. Fig. H shows a plot of the crack tip velocity as 
a function of the energy per unit length (.F/3^)strain in 
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X (in units of w/h) 

FIG. 4: A section of the system showing the (j) = 0.5 level sets 
at equally spaced time intervals. The circle in the middle is 
the initial "crack". After an initial slow transient the crack 
reaches a constant velocity. When it reaches the boundary 
and senses its periodic image, the speed will increase again 
as the two tips coalesce and form a continuous crack in the 
{/-direction (not shown here). This simulation corresponds to 
the data point in Fig. |^ with (.7^/3^)strain ~ 1 and v ~ 0.03. 



the uncracked region. A quadratic fit of [T /y)sf[a\n as 
a function of the crack velocity shows that the fracture 
threshold is (J^/3^)strain = 0.247 ± 0.003. 

A comment should be made about our use of periodic 
boundary conditions. Unhke the case of dynamical frac- 
ture, our overdamped dynamics do not suffer from elastic 
fields reflected from the boundaries or impinging from 
periodic images. However, the periodic boundary con- 
ditions do have important effects. As mentioned above, 
the velocity of the crack tip changes when the periodic 
images of the crack get sufficiently close to each other. In 
addition, for sufficiently thin systems (small width in the 
x-direction) and low strains (w ~ 0), the energy released 
due to the Poisson ratio is enough to favor a crack grow- 
ing in the horizontal direction, parallel to the a;-axis (lim- 
iting how thin we can make our rectangular region) . We 
have minimized the effects due to the periodic boundary 
conditions by using a wide system and by only measuring 
the crack velocity when the tip is far from the bound- 
aries, though further measures might be needed for more 
sophisticated simulations. Note that periodic boundary 
conditions make it hard to apply stresses to the system; 
all the simulations presented in this paper were driven 
by applying an initial strain (equivalent to applying a 
displacement). 
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FIG. 5: The graph shows the crack tip velocity against the 
strain energy per unit length (.F/J^) strain in front of the crack 
tip. The simulations were done with double ended cracks in 
a sample of size X = y — 200, with A = 2/x = 2 and an 



initial strain in the y-direction ey 



0. The extrapolation 



of the velocity brings it to zero at a strain energy of about 
(.7^/3^)strain = 0.247 ± 0.003, as compared to the theoretical 
Griflith's threshold of (JP/3^)surfacc = 2/9 fa 0.222. 



Eq. ( p7| ) gives the analytical surface energy of a sta- 
tionary interface. With tyy = 0, A = 2/i = 2, the total 
free energy per length far behind the crack tip due to 
the two interfaces is[|o| (:?^/3^)surfaco = 2/9 « 0.222. If 
our model were to obey the one proposed by Griffith, 
the crack tip velocity should tend to zero as {T /y) strain 
approaches (.F/3^)surtacc- Apparently, there is a disparity 
between the numerical fracture threshold and the analyt- 
ical Griffith's threshold for our model. Most of the energy 
in front of the crack tip is transferred to the newly cre- 
ated crack surfaces, but some extra energy is needed to 
drive the crack forward. 

Upon investigation, we find a large residual strain Cyy 
building up across the width of the sample near the height 
of the crack tip during fracture. In front of the crack, the 
material has contracted in the y-direction, presumably 
because of the Poisson ratio and the positive strain per- 
pendicular to the crack. Behind the crack, the material 
has been stretched in the y-dircction; this is necessary 
to make up for the compressed material in front of the 
crack. 

It turns out that the difference in energy between the 
fracture threshold and the Griffith's threshold roughly 
matches the residual energy stored in this deformation. 
Typically, viscous effects vanish as the tip velocity goes 
to zero, but we claim that this one goes to a constant 
because the width of the deformation diverges in the limit 
of a stationary crack tip. Below we will give a crude 
analytical argument to show that this is not a finite size 
or finite velocity effect, and give a rough estimate of the 
missing energy. 

Consider an infinitely long system in the y-direction, 
but with a fixed width X . To make the system one dimen- 
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sional, we will ignore any variations in x. We introduce 
here a frame that is stationary with respect to the crack 
tip, y — y ^ vt and t = t, which means that dt — > —vdy 
since any derivatives with respect to t disappear for a 
stationary solution. 

Both far in front of and far behind the tip we assume 



that ey 



0, which means that the displacement field 



Uy{y) is constant in these areas. In the deformed region in 
between (that is, where the crack tip is), the displacement 
field changes by a value Auj^ — Uy(— oo) — Uy(po) (notice 
that the displacement field has a lower value in front of 
the tip). According to Eq. (|^), the strain far ahead of 
the tip exerts a stress 



^yy — XCxx 



(29) 



on the deformed region. Since there is no strain far be- 
hind the tip, this stress must be countered by the viscous 
force due to the movement of the displacement field: 



dUy 



dUy 

dy 



(30) 



Knowing how much Uy changes, we want to find its 
shape in the deformed region. Using that d^Uy = and 
dfUy = —5T/5uy = (A + 2^)dyUy, we have that 



d 
dt 



V 



dUy 

~dt 



dUy 

dy 



du,. 



= {X + 2^i)^ + = 



dy 



dy 



(31) 



A solution to this second order differential equation apart 
from the crack tip position y ~ vt is 



uy{y,t) 



AUy , 

Auy exp 



X+2fj. 



[y - vt) 



y < vt 
y > vt 



(32) 



where we have used that Uy(oo) = 0. 

We can now find the energy dissipated due to the ma- 
terial deformation around the crack tip. Ignoring the 
phase- field in Eq. (|l0|), our crude estimate gives the en- 
ergy dissipated per unit length as the inverse velocity 
times the energy dissipated per unit time: 
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V dt 
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dy 



iXtxxfX 

(A + 2fx)2 ' 



(33) 



where we used Eq. ( |30| ) in the last step. (To quadratic or- 
der in €xx, this equals A^/4/x(A -I- 2fi) times {J^ /y)sti-aim 
so the dissipation is a fixed fraction of the total avail- 
able strain energy. For our system, the fraction is 1/8.) 
We thus conclude that the energy dissipated due to the 
deformation around the crack tip only depends on the 
initial strain and the width of the system. In our sys- 
tem of width X — 200, the strain corresponding to the 
fracture threshold is about exx = 0.026. For numerical 
comparison, we insert the strain and the other values 
that we used in our simulations into Eq. (p3|), and find 
that T/ydiss ~ 0.07, which is about twice the difference 
between J^/3^strain and J"/3^surtaco • Even though we ig- 
nored both the phase-field and all variations in x, our 
crude estimate shows that as the velocity goes to zero a 
fixed fraction of the energy in front of the crack is used 
to strain the material in a diverging region around the 
crack tip before it is dissipated. 

We also did a simulation to verify that Eq. ( ^ ) is of 
the right form. With {X,y) = (100, 1200) and Cxx = 0.08 
(which gives (.F/3^)strain — 1-01), we did see an exponen- 
tial decay of Uy in front of the crack tip. The numerical 
values were slightly off, with the exponential decay being 
twice as fast as predicted. This was to be expected since 
we have ignored the phase-field completely, which clearly 
couples to the displacement field at long distances. Our 
simulation area still was not long enough for Auy to sat- 
urate, but it reached about 2/3 of the value predicted in 
Eq. @). 

In Fig. H, we plot the velocity as a function of strain 
energy available for fracture (.F/3^) strain- In the fracture 
literature, the mode I fracture threshold is usually quoted 
in terms of the energy release ratepT]: the strain energy 
Q flowing into the crack tip. This can be measured, for 
example, by the use of J-integrals|jl^, see Appendix |C[ 

From the discussion above, we conclude that the dis- 
parity between the Griffith's prediction and the fracture 
threshold as measured by (.F/3^)strain is due to energy 
dissipation far from the crack. Thus we cannot expect 
G to equal (.F/3^) strain- Indeed, the fracture threshold 
measured using a J-integral close to the crack tip should 
agree with Griffith's prediction. 



V. FUTURE WORK 

Fracture often occurs in crystalline materials, and it 
therefore seems reasonable to add anisotropy to quanti- 
ties such as mobility, surface tension, and elastic strain 
energy. There might also be unwanted anisotropy in- 
cluded due to the underlying grid in the numerical 
solver, and one might be able to repress this effect by 
adding "counter" -anisotropy in the quantities mentioned 
above. In addition to anisotropy, one could add some 
noise; either spatially in the Hamiltonian or initial con- 
ditions, which would break the symmetry and add hetero- 
geneities, or temporally, which would mimic the effects 
of fluctuations due to a finite temperature. 
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Because of the fourth-order gradients, our simulation 
is rather numericaUy intensive. Using a multiresolution 
grid would speed up the calculation Q. The main idea 
is that additional information can be added between ex- 
isting grid points where a linear interpolation would give 
poor results (that is, add grid points locally where the 
solution has high frequency components). This way one 
can achieve a solution with uniform accuracy and avoid 
unnecessary use of computational resources. A similar 
approach can be found in[p^, where an adaptive mesh in 
a phase-field solidification problem is used to add detail 
only at the boundaries where it is needed. 

Currently, the numerical integration of the phase-field 
equations is done partly using Fourier transforms. This 
allows us to do all the linear parts of the equations implic- 
itly without solving huge linear systems. On the other 
hand, the current implementations puts several limita- 
tions on the kind of simulations that can be done; es- 
pecially, it is hard to do simulations with other than 
periodic boundary conditions. It thus seems beneficial 
to use real-space implicit methods (and perhaps add in 
multi-resolution capabilities at the same time). 

This paper is only concerned with mode I fracture in 
two dimensions. It would seem reasonable to try to ex- 
tend the model to include mode II and mode III as well. 
Mode II could be done either by shearing the model with 
a crack running vertical or horizontal, or by squeezing 
one way and extending the other with a strain s and 45 
degree crack, where s is equal to the shear strain divided 
by v^. One way to start exploring mode II could be to 
use finite element calculations to get an initial sheared 
configuration, and then use the phase-field model to re- 
lax this system. Mode III has recently been explored 
in Ref. Q, where a different, non-conserved phase-field 
model is used. Ultimately, the goal is to do three di- 
mensional simulations incorporating all three modes of 
fracture, given the success of the two-dimensional mod- 
els. 
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APPENDIX A: REDUCED UNITS 

Most of the equations in this paper are written in unit- 
less form, where the basic units are: 



I = w/h 
d = 1/r] 



(length) (Ala) 
(energy density) (Alb) 
(diffusivity) (Ale) 



Here w'^ is the cost of gradients of (f> in Eq. (|l|), h^/6A 
is the height of the energy barrier between the phases in 
Eq. (H), and 77 is the viscosity controlling the response of 
the displacement field in Eq. (8b). The equations were 
made unitless by replacing all quantities (say, x and t) 
by their unitless counterparts multiplied by the appro- 
priate combination of basic units (as in x*l and t*P / fd, 
* = unitless), and then choosing the basic units Z, / and 
d conveniently. To transform the quantities back from 
unitless form, just multiply them with their basic unit 
{x = x*l). 

As a reminder, the variational derivatives of the free 
energy and the Lame constants A and /i all have units 
of energy density (/), and the diffusion constant D has 
units of diffusivity (d). 



APPENDIX B: LAME CONSTANTS 

We have to check what the relation between A, /i and 
v is in our model: the addition of (f> changes the effective 
elastic constants on long wavelengths from the values in- 
put to the free energy in Eq. (^. In the following, two 
and three dimensional refer to the mathematical dimen- 
sionality of the system, not special cases like plain strain. 
Start with a rectangular sheet (2D) or block (3D) of ma- 
terial, with height h in the y-direction, and width w in 
the a;-direction (and in the z-direction if it is in 3D). The 
idea is then to strain the material an infinitesimal amount 
in the y-direction until the height becomes /i', and min- 
imize the free energy with respect to the new width w' . 
The Poisson ratio v is then given by 



lim 

h'^h 



-yy 



(Bl) 



Before doing the model proper, we start off with a free 
energy that only includes the elastic energy: 



= / dV£[e 
Jv 



(B2a) 



In the scenario that is described above there is no shear, 
so eij = if i ^ j. Further, eyy = {h' — h)/h', and 
^xx — ^zz = {w' — w) / w' (in 2D, Czz is non-existent). The 
integrand is constant everywhere, so the integral turns 
into a fa ctor h'jw')'^ (or h'w' in 2D). This is inserted 
into Eq. (B2a). Next, one finds the minimum of J-q with 
respect to w' (keeping all the other variables constant). 
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Dimension 






To 








2 


hw 
h' w' 






A = 


2^1^ 1 
1-1/ 2 


1 — 1/ 


2 


1 


A 


2^1/ 

^ 1-v 


A 


2^1/ 

^ 1-1/ 


\ _ 2^1^ + 1 1 
1 — 1/ 2 


3 


1 


A 


2^1/ 

^ l-2i/ 


A 


2fj,h' 

^ l-2i/ 





TABLE I: Relations between A, ^ and v when a material is 
stretched infinitesimally. 



and inserts into Eq. (Bl), taking the limit. Solving for A 
does indeed give the standard answers, for both two and 
three dimensions. 

The next two energies that were tried, were 



•^1 = / dv\. 



If + [e] , 



(B2b) 
(B2c) 



The gradient term is zero, since (/> is uniform. Having 
introduced </), we need to add another restriction; con- 
servation of mass is given by (pohw'^ = (f)'h'{w'f (or 
4>oh'w = (p'h'w' in 2D). Here, (j)o is the value when the 
material is relaxed, while (j)' is the value after the ma- 
terial has been stre tched [it is the latter which will be 
inserted into Eqs. ( B2b| ) and (B2c)]. We assume that 
(j>Q = 1. We checked two cases for 0'; one where (/>' = 1 
(the density was not allowed to change), and one where 
4>' = hw/h'w' (the latter is the 2D expression; having (p 
change in three dimensions was quite hard to calculate). 

The results after minimizing the free energy, finding 
taking the limit and solving for A can be found in Table |. 
As one can see, adding the phase-field cf) to the model also 
requires that the double-well potential includes the e,„m 
term (through (/'^[e]) in order to get the right relation 
between the material constants. 

A final comment: The model described in this paper is 
two-dimensional, so the maximum range for the Poisson 
ratio in two dimensions is Q < v < 1. This is equivalent 
to a three dimensional model under plain strain, where 
the corresponding range for the Poisson ratio is < < 
1/2. Thus the latter expression (for three dimensional 
plain strain) is used throughout the paper, often using 
the value v — l/i. 



APPENDIX C: THE J-INTEGRAL 

The energy release rate Q can be calculated for our 
problem (where the crack is parallel to the y-axis) using 
the Jy component of the J-integral. Instead of perform- 
ing the line-integral, it is common to convert it to an 
area-integral for increased accuracy when doing the inte- 
gral numerically. The area integral is defined as 



Jy 



Q.{x, y)dxdy 



where 



dy dxj 



(CI) 



(C2) 



Here q is a function that is unity around the crack tip 
and zero outside. Notice that if q is constant in a re- 
gion, n{x, y) = 0, so in effect the line integral is replaced 
by a "thick line" integral, where the "thick line" exists 
everywhere q has a gradient. 

A small complication is that our f2(x, y) is given in de- 
formed coordinates £,x{x, y) = x + u^ix, y) and £,y{x, y) = 
y + Uy{x,y): 



Jy = ^ m.,Cy) \Ji^.,Lj)\ d^xd^y (C3) 
J A 

where is the Jacobian given by 



dx dx 
dy dy 

J 



(C4) 



Using the identity J{^x,Cy) ■ J{x,y) = /, one gets 
1 



1 



dy Oy ^y 

1 

(1 + £xx){l + £yy) - dyUxdxUy 



(C5) 



The ener gy release rate can thus easily be computed using 
Eq. (C3), where S^x and are the usual coordinates on 
the phase-field grid. 
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If we had chosen a non-zero constant of integration in 

Eq. (0), then Eq. (^) would not have been linear in 

(j), and T[(j>] would not have been a simple fourth order 

polynomial. 

For easy reference, this means that A -f 2^ = -BV(1 - 

2B^), and if A = 2^i = 2 then B = ±2/3. 

(j) = is also a local minimum of the free energy density, 

but this corresponds to the gas phase. 

After relaxing the interface numerically, the actual value 

of the surface energy is {T/y) surface 

= 0.219 ± 10"^ 

The reason it is not 2/9 is that the numerical solution 
is slightly different from the analytical one due to the 
discrete nature of the simulation. 

For plane strain, the fracture threshold is given in terms 
of the stress intensity factor Ki, which is related to the 
energy release rate by C/ = K^{1 ~ '^)/2/i. 



